NASA Contractor Report 187629 
ICASE Report No. 91-75 


ICASE 


AN APPROXIMATE RIEMANN SOLVER FOR 
HY PERVELOCITY FLO WS 


(NASA— CR— 187629 ) AN APPROxlHAT F R I F MANN N91-32463 

SILVER FOR HY n ERVFLOCTTY FLOWS Final Report 
(ICASE) 1 6 p CSCL 20D 

Unci as 

G3/34 00A67B8 


P. A. Jacobs 


Contract No. NAS 1-1 8605 
September 1991 

Institute for Computer Applications in Science and Engineering 
NASA Langley Research Center 
Hampton, Virginia 23665-5225 

Operated by the Universities Space Research Association 


NASA 

National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton, Virginia 23665-5225 



AN APPROXIMATE RIEMANN SOLVER 
FOR HYPERVELOCITY FLOWS 


P. A. Jacobs 1 

Institute for Computer Applications in Science and Engineering 
NASA Langley Research Center 
Hampton, VA 23665 


ABSTRACT 

We describe an approximate Riemann solver for the computation of hypervelocity flows 
in which there are strong shocks and viscous interactions. The scheme has three stages, 
the first of which computes the intermediate states assuming isentropic waves. A second 
stage, based on the strong shock relations, may then be invoked if the pressure jump across 
either wave is large. The third stage interpolates the interface state from the two initial 
states and the intermediate states. The solver is used as part of a finite-volume code and is 
demonstrated on two test cases. The first is a high Mach number flow over a sphere while 
the second is a flow over a slender cone with an adiabatic boundary layer. In both cases the 
solver performs well. 


1 Research was supported by the National Aeronautics and Space Administration under NASA Contract 
No. NAS 1-1 8605 while the author was in residence at the Institute for Computer Applications in Science 
and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23665. 
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Nomenclature, Units 

a : local speed of sound, m/s 

E : total energy (internal + kinetic), J/kg 

e : specific internal energy, Jjkg 

h : specific enthalpy, J/kg 

M : Mach number 

P : pressure, Pa 

Pr : Prandtl number, ( C p p/k ) 

R : gas constant, J / kg j K 

Re : Reynolds number 

T : temperature, I( 

t : time, s 

U : Riemann invariant 

u : ar-component of velocity, m/s 

v : y-component of velocity, m/s 

ws : wave speed used in the Riemann solver 

x : x-coordinate, m 

y : y-coordinate, m 

Z : intermediate variable 

a : weighting function 

p : density, kg/m 3 

7 : ratio of specific heats 

p. : coefficient of viscosity, Pa.s 

Superscripts 

* : intermediate states for the Riemann solver 

/locally tangent to the cone surface 

Subscripts 

MIN : minimum allowable value 

e : boundary-layer edge condition 

x,y : cartesian components 

L y R : left state, right state 
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1 Introduction 


In recent years the proliferation of relatively fast computers has popularized the direct cal- 
culation of viscous, compressible flows in a time-accurate manner. In some situations, such 
as the transient hypersonic flow over a model in a shock-tunnel, numerical simulation is the 
only way to extract detailed information about the flow field. Such computations are very 
demanding as there are both strong shocks and rarefactions and strong viscous interactions. 

This note describes a robust Riemann solver for use in transient hypervelocity flow calcu- 
lations. The full code [l] is based on a cell-centred time-dependent finite- volume formulation 
of the axisymmetric Navier-Stokes equations in which the governing equations are expressed 
in integral form over arbitrary quadrilateral cells. The time rate of change of conserved 
quantities in each cell is specified as a summation of the fluxes through the cell interfaces. 
The inviscid components of the fluxes are computed with the approximate Riemann solver 
while the viscous fluxes are calculated by application of the divergence theorem. At each 
time step, we first interpolate the flow state (consisting of a set of values for p , u, w, e, P , a) 
to either side of each interface at the start of the time step and then apply a Riemann solver 
to estimate the flow state at the interface during the time step. Note that the solver is 
applied in a locally rotated frame of reference in which the u-velocity is normal to the cell 
interface. 

There are a number of Riemann solvers that can be used including exact iterative 
schemes [2] and approximate (noniterative) schemes [3, 4]. The approximate schemes are 
generally less computationally expensive than the iterative schemes and, because the Rie- 
mann solver comsumes a large fraction of the total computational effort, an approximate 
scheme is favoured. Although the Roe-type solver is popular because it is relatively fast, 
there are situations such as the double- Mach- reflection case [5] and flow over a sphere [6] 
where it may produce spurious results. The Osher-type solver [4] is considered to be fairly 
robust and free of adjustable parameters, however, we have experienced some difficulty in 
applying it to flows with very strong shocks. 

Here, we take a middle road between the fully iterative solvers and the single-step ap- 
proximate solver and, at the cost of some computational expense, produce an approximate 
solver which is reliable in extreme flow situations and is vectorizable with current compilers 
for vector computers. 
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2 Approximate Riemann Solver 


The current solver is a 3-stage approximate Riemann solver in which the first stage computes 
the intermediate pressure and velocity assuming isentropic wave interaction. A second stage, 
based on the strong-shock relations, may be invoked to improve the first-stage estimate if 
the pressure jump across either wave are sufficiently large. In practice, this modification 
has been required only in extreme conditions such as those found in the bluff-body test case 
(Section 3.1). The final stage is to select /interpolate the interface state (p, u, v, e, P, etc) 
from the set of left, right and intermediate states. If stage 2 (strong shock modification) is 
not invoked, the solver is much like Osher’s approximate Riemann solver [4]. 


STAGE 1: The first stage of the Riemann solver assumes that a spatially constant left state 
(subscript L) and right state (subscript R) interact through a pair of finite-amplitude (and 
isentropic) compression or rarefaction waves. Perfect gas relations ([7] cited in [2]) are used 
to obtain the intermediate states (L*, R*) in the gas after the passage of left-moving and 
right-moving waves, respectively. The expressions implemented in the code are 


P*. = Pr = P* = Pi 
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where the Riemann invariants are 
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and the intermediate variable Z is given by 


= an /P^G-m^) 
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Note that these expressions involve the power operator which is computationally expensive. 
For a limited range of base and exponent, the standard power function is replaced by the 
approximate expansion [1]. In the exceptional situation of (Ul — Ur) <0, we assume that 
a (near) vacuum has formed at the cell interface and set all of the interface quantities to 
minimum values. 


STAGE 2: If the pressure jump across either wave is large (say, a factor of 10), then the 

guess for the intermediate pressure is modified using the strong shock relations. 
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If P* > 10 Pl and P* > 10 Pr then both waves are taken to be strong shock waves and 
the intermediate pressure and velocity can be determined directly as 


p* t "f” i _ 

1 - o Pl 


\/PR 
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{u L - Ur) 
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and 
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( 6 ) 


If P* is greater than Pi or Pr (but not both), the stage-1 estimate for P* can be improved 
with two Newton-Raphson steps of the form 


where 


and 
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( 10 ) 


The strong-shock expressions used in (9), (10) can be obtained from the normal shock ex- 
pressions in [8] by taking the limits as the pressure jump becomes large. During the update, 
we ensure that P* > Pmin where Pmin is some small value. After updating P*, the inter- 
mediate velocity is evaluated using the relevant strong-shock relation from (9) or (10). 


STAGE 3: Now that we have computed the pressure and velocity in the intermediate regions 
behind the waves, the other intermediate flow properties may be evaluated. The interface 
conditions used in the inviscid flux vector can then be selected or interpolated from the 4 
flow states using the logic shown in Fig. 1. Note that although only the left-moving wave is 
discussed below, a similar procedure is used to obtain the flow state behind the right-moving 
wave. 


If the pressure rises across the left-moving wave (i.e. P* > Pi ), the left wave is assumed 
to be a shock and density is obtained from the Rankine-Hugoniot relation as 


(7 + 1)-P* + (7 ~ 1)Pl 
X l + 1 )Pl + {l ~ 1)P\ 


( 11 ) 
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The specific internal energy is obtained from the equation of state as 

P* 


( 12 ) 


‘ L (7 - IK ’ 

and estimates for the local speed of sound (for later use in the interpolation of the interface 

(13) 


properties) are 


(14) 


«l = \AKt - ije£ • 

The velocity of the wave (relative to the initial left state) is given by 

h+IPL (P- . 7 -lM ‘ /2 

u L - ws L = — — + — I 

L 2 p L \P L 7 + 1/J 

where wsl is the velocity of the wave relative to the cell boundaries. 

If the pressure falls across the left-moving wave (i.e. P * < Pl), the isentropic-wave 
relations are used to obtain the intermediate properties. The local speed of sound is obtained 
from the Riemann invariant as 


al = (U L -ul)( 7 -l)/2 , 

while the specific internal energy is obtained from the sound-speed relation as 

, K ) 2 

(7-1)7' 

The density is obtained from the equation of state as 

. P * 

Pl - 


(15) 


(16) 


(17) 


(7 - 1 K 

and the velocity of the leading-edge of the wave (relative to the initial left state) is given by 

(18) 


Ul — WSL = a L 


3 Test Cases 

3.1 High Mach Number Flow around a Sphere. 

The robustness of the code is demonstrated by computing a Mach 60 flow over a 1 .hmm 
radius sphere with a domain consisting of a 60 x 60 mesh of cells. The y = 0 boundary. is 
the symmetry line (and stagnation line) while a tangency condition is applied at the surface 
of the sphere. Free-stream conditions of 

p = 0.5097 x 10" 2 kg/m 3 , P = 427.1 Pa, e = 2.095 x 10 5 J/kg, 
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u = 20600 m/s , v = 0, M nom inai = 60.1 


are applied to the curved inflow boundary, the shape of which is derived from the shock- 
position correlations in [9]. Flow conditions at the outflow boundary are obtained by zero- 
order (constant) extrapolation. Initial conditions throughout the domain are set to 

p = 0.5097 x 1(T 2 kg/m 3 , P = 427.1 Pa, e = 2.095 x 10 5 J/kg, u = 0, v = 0. 
Despite the very high temperatures in the shock layer, the gets is considered perfect with 

7 =1.4, R = 287 J/kg/K, Pr = 0.72, 

and Sutherland’s law is used for the coefficient of viscosity. The Navier-Stokes equations 
are then integrated forward in time using high-order MUSCL interpolation and Euler time- 
stepping with a CFL number of 0.5. 

Figure 2 shows the flow field (pressure and Mach contours) at t = 13.6/xs after the flow 
has approached steady state. Discrete points from experimentally derived correlations [9] 
are plotted on the pressure contours. Given that M — 60 is beyond the range of the data 
used for the correlations, agreement is good. The largest deviations are near the outflow 
boundary. Profiles of density and pressure along the line of cells adjacent to the x-axis are 
shown in Fig. 3. The shock appears to be captured in 2 or 3 cells with no oscillation and the 
density jump is close to the ideal strong-shock value of 6. The pressure ratio from free-stream 
to the stagnation point is 4621 which is very close to the ideal value of 4636 for M — 60 (see 
e.g. [8] Table II). 

A similar calculation with M = 12 was reported in [1] and, for that condition, an Osher- 
type solver (i.e. stages 1 and 3 only) failed to produce a solution. Also, a finite-difference 
scheme using Roe-type flux-difference splitting required a rather large value for its entropy- 
fix parameter in order to obtain a physically reasonable solution (J. White, NASA Langley 
Research Centre, private communication). Also note that, while the M = 60 shock is very 
strong, the high temperature in the region behind the shock enhances the viscous dissipation 
and may result in a smoother solution than seen at lower Mach numbers. 


3.2 Flow Over a Sharp Cone 

To illustrate the behaviour of the solver in the presence of strong viscous effects, we show 
the computed results for M ~ 8 flow over a sharp 7° cone. The axis of the cone is aligned 
with the free stream. 
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Two cases are considered in which the cone flow domain is discretized as a set of 100 x 60 
cells and 100 x 90 cells. Free-stream conditions of 

p = 1.0809 x 10“ 2 kg/m 3 , P = 165.51 Pa, e = 3.8281 x 10 4 J/kg, 

T = 53.35 I<, u = 1164.0 m/s, v = 0, M nomina , = 7.95 

are applied to the left and upper boundaries while the outflow boundary conditions are 
obtained by extrapolation and the cone surface is modelled as a no-slip, adiabatic boundary. 
To match the experimental conditions in [10], the gas was considered to be a perfect gas 
with 

7 = 1.4, 77 = 287 J/kg/K, Pr = 0.7, 
and viscosity was obtained from the Sutherland expression 

7 * 3/2 

" = 1611 * 10 rTTuS p “ s ' 

Based on free-stream conditions and the length of the cone, the Reynolds is approximately 
3.3 x 10 6 . The initial state of the flow in the domain is 

p = 1.0809 x 10~ 2 kg/m 3 , P = 165.51 Pa, e = 3.8281 x 10 4 J/kg, u = 0, v = 0. 

The Navier-Stokes equations are then integrated forward in time using high-order MUSCL 
interpolation and Euler time-stepping with a C FL number of 0.5. 

Figure 4 shows the flow field (pressure and density contours) t = 22 ms after the flow 
has approached steady state. The pressure field is almost conically symmetric, as per the 
inviscid solution of Taylor and Macoll [11] (see also [12] Ch. 10) and the shock angle is 
still approximately the inviscid value of 10.5° ([8], Chart 4). The shock, however, is slightly 
curved near the apex of the cone. Boundary layer profiles of velocity and temperature at 
x* = 1.0m are shown in Fig. 5 for both the present finite- volume solutions and a boundary- 
layer solution using edge conditions of 

p e = 2.044 x 10~ 2 kg/m 3 , P e = 416.7 Pa, u* = 1148.6 m/s, T e = 71.04 K. 

There is good agreement between the present finite-volume solutions and the spectrally- 
accurate solution [13], especially near the cone surface. Although the outer region of the 
boundary layer is underresolved, (even for the 100 x 90 mesh) the finite-volume solutions 
appear to be converging to the spectral solution. 
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if («* > 0) then 


The contact-discontinuity has moved to the right 
and the interface state is determined from the 
L and L* states. 

if (P* > Pl) then 
The left-moving wave is a shock, 
if ( wsl > 0) then 

All waves have moved to the right. 

Interface state is equal to L. 
else 

Interface state is equal to L*. 
endif 
else 

The left-moving wave is a rarefaction, 
if («l - ol > 0) then 

All waves have moved to the right. 

Interface state is equal to L. 
elseif (u* L — a* L > 0) then 
The rarefaction straddles the interface. 
Interpolate the interface state from 
states L and L*. 
else 

The entire rarefaction moved to the 
left of the interface. 

Interface state is equal to L*. 
endif 
endif 

else 

The contact discontinuity has moved to the left 
and the interface state is determined from the 
R and R* states in a similar manner... 


Figure 1: Interpolation logic for the Riemann solver. 
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Figure 5: Comparison of the present finite- volume solution with a spectral solution at 1.0m 
from the cone apex: (a) tangential velocity; (b) temperature. Solid line represents spectral 
solution, (3 denotes 100 x 60 mesh, A denotes 100 x 90 mesh. 
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